Mon. Not. R. Astron. Soc. 000. [TUTUl (20051 Printed 2 February 2008 (MN WF&L style file v2.2) 

Cusps in CDM halos 



Jiirg Diemand 1,3 *, Marcel Zemp 1 ' 2 , Ben Moore 1 , Joachim Stadel 1 & Marcella Carollo 2 

1 Institute for Theoretical Physics, University of Zurich, Winterthurerstrasse 190, CH-8057 Zurich, Switzerland 

2 Institute of Astronomy, ETH Zurich, ETH Honggerberg HPF D6, CH-8093 Zurich, Switzerland 

3 Department of Astronomy and Astrophysics, University of California, 1156 High Street, Santa Cruz CA 95064, USA 



2 February 2008 



in 
o 
o 

(N 

Oh 

m 

(N 

m 
> 
in 

(N 
O 

in 
o 

Oh 

6 

c3 



13 



ABSTRACT 

We resolve the inner region of a massive cluster forming in a cosmological ACDM 
simulation with a mass resolution of 2 x 10 6 Mq and before z=4.4 even 3 x 10 5 M Q . 
This is a billion times less than the clusters final virial mass and a substantial increase 
over current ACDM simulations. We achieve this resolution using a new multi-mass 
refinement procedure and are now able to probe a dark matter halo density profile 
down to 0.1 percent of the virial radius. The inner density profile of this cluster halo 
is well fitted by a power-law p oc r~ 7 down to the smallest resolved scale. An inner 
region with roughly constant logarithmic slope is now resolved, which suggests that 
cuspy profiles describe the inner profile better than recently proposed profiles with 
a core. The cl uster studied here is on e out of a sample of six high resolution cluster 
simulations of iDiemand et all l)2004bj) and its inner slope of about 7 = 1.2 lies close 
to the sample average. 

Key words: methods: N-body simulations - methods: numerical - dark matter - 
galaxies: haloes — galaxies: clusters: general 



1 INTRODUCTION 

Recently a great deal of effort has gone into high reso- 
lution simulations which have revealed density profiles of 
cold dark matter halos down to scales well below one per- 
cent of the virial radiu s jFu4cushige J _i<awaQ^ 
i Tasitsiomi et all20o4lNavarro et all2004tlReed et all2005t 
biemand. Moore fc Stadelll2004bl "DMS04" hereafter). But 
the form of profile below ~ 0.5 percent of the virial radius 
remained unclear and there was no clear evidence for a cusp 
in the center, i.e. no significant inner region with a constant 
logarithmic slope. Galaxy cluster halos would be the ideal 
systems to resolve cusps numerically because of their low 
concentration. In a galaxy or dwarf halo the inner power 
law is much harder to resolve because it lies at a smaller 
radius relative to the size of the system. 

The existence of a core or a cusp in the center of 
CDM halos has important observational consequences and 
is the crucial point in many tests of the CDM theory. Com- 
parisons of dark matter simulations to rotation curves of 
low surface brightness galaxies (LSB) seem to fav or con- 
stant density cores for most ob se rved systems fe.g. iMoore l 
20041: iFlores fc Primackl Il994l : ISalucci fc Burkerd l2000t 
deBl ock et Id! 2001; see, however Ivan den B osch fc Swatersl 
20011; IS waters et alJ Eoolt ISimon et alJ l2005h x. But these 
comparisons still depend to some extend on extrapolations 
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of the simulated profiles toward the center: IStoehrl d2004h 
extrapolate to a constant density core and claim that the 
discrepancy to LSB galaxy rotation curves is much smaller 
than previously believed. 

The strength of the 7-ray signal from dark matter 
annihilation depends on the square of the dark matter 
density and the calculated flux values spread over sev- 
eral orders of magnitude, depending on how one extrapo- 
lates the density profiles from the known, resolved regions 
down into the centers of the galactic halo and its sub- 
halos JCalcaneo-Roldan fc Moorell2000[ IStoehr et alJl2002l : 



iBertone fc Merritl 120051 iPrada et alJ l2004h . Small, very 
abundant, Earth to Solar mass su bhalos could be very lu- 
minous in 7-rays if they are cuspy JPiemand et alJl2005h . 

The highest resolutions in cosmological simulations are 
reached with the widely used refinement procedure (e.g. 
lBertschinge'ill200lD : First one runs a simulation at uni- 
form, low resolution and selects halos for re- simulation. 
Then one generates a new set of initial conditions using 
the same large scale fluctuations and higher resolution and 
additional small scal e fluctuations in the selected r egion. 
With this technique iNavarro. Frenk fc White) ^99(S) were 
able to resolve many halos with a few ten thousand parti- 
cles and to infer their average density profile which asymp- 
totes to an p(r) oc r _1 cusp. Other authors used fittin, 
funct i ons with steeper ( - 1.5) cusps j Fukushigc & Makin 
1997t iMoore et all Il998t iMoore et all Il999t iGhigna et al 



2000). Small mass CDM halos have higher concentra- 
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tions due to their earlier collapse iNavarro et all Il996h 
but the slopes of the i nner density profil e s are inde- 
pend ent of halo mass ijMoore et alJ 1200 it ICoh'n et alJ 
2004). Open, "standard" and lambda CDM cosmologies, 
i.e. models with (fi M) fi A ) = (0-3 0.0), (1.0,0.0) and 
(0.3, .7) yield equal inner profiles jFukushige fc Makinol 
l2003l: iFukushiee. Kawai fc Makinol l2004h . There is some 
indication that models with less small scale power like 
WDM lead to sh allower inner profiles (e.g. lCoh'n et all200ft 
iReed et all2005l) . Different equation of states of the dark en- 
ergy component lead to different collapse times and halo con- 
centrations but it is not cl ear yet if it also affects slopes well 
insid e of the scale radius jMaccio et aljl2004t iKuhlen et alJ 
120041) . Most current simulations do not resolve a large 
enough radial range to determine both the concentration 
and the inner slope; at the c urrent resolution th ese param- 
eters show some degeneracy jKlvpin et alJl200lj) . 

Recently a large sample of ACDM halos resolved 
with a million and more particles was simulated 
JSprineel et al.ll200ldlTasitsiomi et aJl2004 iNavarro et alJ 
120041: IReed et alJ 120051: iGao et alJ I2005T) and the best 
resolved systems contain up to 25 million particles 
JFukushige. Kawai fc Makinol 12004 DMS04). But even 
these very large, computationally expensive simulations re- 
solved no inner region with a constant l ogarithmic sl ope. 
jNavarro et alJ 120041 : IStoehr et ail 120021 IStoehJ I2004T) in- 
troduced cored profiles which seem to fit the simulation 
data better than th e cu spy profiles propos ed earlier by 
iNavarro etafl Jl996ft and lMoore et alJ ll99Sft . This better 
fit was interpreted as indication against cuspy inner pro- 
files. However these cored profiles have one additional pa- 
rameter and therefore it is not surprising that they fit the 
data better. DMS04 showed that an NFW-like profile with 
the inner slope as additional free parameter fits the highest 
resolution profiles just as well as cored profil es. Some the- 
orctical argument s seem to favor cusps (e.g. iBinnevlliool 
lHansen fc Moorel l2004h but make only vague predictions 
about the inner slopes. A recent model combines simula- 
tion results a nd analytical argume nts to predict an inner 
slope of -1.27 dAhn fc Shapiroll2005l) At the moment higher 
resolution simulations seem to be the only way to decide the 
core vs. cusp controversy. 

Here we present simulations of one of the galaxy clus- 
ters from DMS04 with two orders of magnitude better mass 
resolution. Our results give strong support to cuspy inner 
profiles. This increase in resolution was made possible with 
only a moderate increase in computational cost by using a 
new multi-mass refinement technique described in Section 
121 In Section |3 we present our results and in Section 0] the 
conclusions. 



2 NUMERICAL EXPERIMENTS 

Table gives an overview of the simulations we present in 
this paper. All runs discussed in this paper model the same 
ACDM cluster labeled "D" in DMS04. With a mass reso- 
lution corresponding to 1.3 x 10 s and 1.04 x 10 9 particles 
inside the virial radius of a cluster, DM25 and DM50 are 
the highest resolution ACDM simulation performed so far. 
Due to the large number of particles and the correspond- 
ing high force and time resolution these runs take a large 



amount of CPU time. Fortunately the inner profiles of CDM 
clusters are already in place around redshift one and evolve 
little between z = 4 and z = JFukushige et al.1 120041 : 
iTasitsiomi et alJl2004HReed et alJl2005h . Therefore one does 
not have to run the simulations to z — to gain insight into 
the inner density profile. We stop DM50 at z = 4.4, DM25 
at z = 0.8 and use the medium resolution runs D5 and D12 
to quantify the low redshift evolution of the density pro- 
file of the same cluster. Run DM25 was completed in about 
2 x 10 5 CPU hours on the zBox supercomputer . The con- 
vergence radius of run DM50 is 1.7 kpc, estimated using the 
r ex A" -1 ' 3 scaling and the measured converged scales from 
DMS04. 

2.1 Multi mass refinements 

Often in cosmological N-body simulations one uses high 
resolution particles only where one halo forms and heav- 
ier particles in the surroundings to account for the external 
tidal forces. One usually tries to defines a large enough high 
resolution region to minimize or avoid mixing of different 
m ass particles with i n the region of interest. One exception 
is lBinnev fc Knebd i2002T) who used particles of two differ- 
ent masses everywhere to estimate the amount of two body 
relaxation in cosmological simulations. In plasma simula- 
tions on the other hand multi mass simulations h ave been 
successfully used since the 1970s (e.g. lDawsonll983l and refs. 
therein). Here we apply this idea to increase the resolution 
in the core of one cluster halo in a cosmological N-body 
simulation. 

The refinement procedure is usually applied to entire 
virialised systems, i.e. one marks all particles inside the virial 
radius of the selected halo and traces them back to the initial 
conditions. Then one refines the region that encloses the po- 
sitions of the marked particles. Usually the region is further 
increased to prevent any mixing of low resolution particles 
into the virial radius of the final system. In DMS04 all par- 
ticles within 4 comoving Mpc in the initial conditions were 
added to the high resolution region. This assures that only 
light particles end up within the virial radius of the final 
cluster and it also has the advantage that halos in the out- 
skirts of the cluster (out to 2 or 3 virial radii) are still well 
resolved jMoore et alJl2004tl . But with this procedure only 
between one fourth to one third of all the high resolution 
particles end up in the cluster. 

If one is only interested in the inner regions of a halo 
it is possible to use a new, more efficient way of refinement: 
Instead of refining the whole virialised system we only refine 
the region where the inner particles come from. This allows 
to reduce the size of the high resolution region considerable, 
because most of particles that end up near the center of the 
system start in a very small region, compared to the region 
which one finds by tracing back all the particles inside the 
virial radius. Using this technique it is possible to reduce 
the computational cost of a CDM cluster simulation by at 
least one order of magnitude at equal force and mass reso- 
lution in the inner region. Of course now one has different 
mass particles inside the final virialised structure therefore 
we must verify that significant equipartition and relaxation 

1 http:/ /www-theorie.physik.unizh.ch/'~stadel/zBox/ 
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Table 1. Parameters of the simulated cluster. At z=0 the viral mass is 3.1 X 10 14 Mq and the virial radius is 
1.75 Mpc. -/Vjjr is the number of high resolution particles and mjjR is the mass and ejjR the force softening 
length of these particles. For the multi-mass runs we also give the masses (ttilr) and softenings (elr) of the 
next heavier particle species. Softening lengths are given at z=0, "[c]" indicates that a constant softening 
in comoving coordinates was used, "[p]" indicates that the softening was constant in physical units after 
z=9 and constant at ten times this value in comoving units before z=9. The resolved scales are constant in 
physical units and give the innermost radius we expect to resolve with the given mass resolution. N v ^ e ff 
is the actual number of particles within the virial radius at z=0 for runs D6, D9 and D12. For the multi- 
mass runs it is the number needed to reach the same resolution in the inner part by doing a conventional 
refinement of the entire system. All runs are 300 Mpc cubes with periodic boundaries, well outside of the 
cluster forming region the resolution is decreased (as in DMS04). 



Run Zstart Zend «HR ^HH WHR e LR "ilR '"resolved V time- N vi[ ^ S 

[kpc] [M ] [kpc] [M Q ] [kpc] step 



DB 


52.4 





4.2[p] 


4'898'500 


D6 


36.13 





3.6[p] 


31'922'181 


DMGse 


36.13 





3.6[p] 


922'968 


DMQle 


36.13 





3.6[p] 


922'968 


D9 


40.27 





2.4[p] 


31'922T81 


DM9 


40.27 





2.4[p] 


3'115'017 


D12 


43.31 





1.8[p] 


14'066'458 


DM25 


52.4 


0.8 


0.84[c] 


65'984'375 


DM25U 


52.4 


0.8 


0.84 [p] 


65'984'375 


DM50 


59.3 


4.4 


0.36[c] 


16'125'000 



jBinnev fc Knebdl2002llDiemand et alJl2004af) is not occur- 
ring and affecting the final results, fn section \l. 31 we show 
that the density profiles of such multi-mass clusters (runs 
DM61e and DM9) are the same as the ones of fully refined 
clusters at equal peak resolution (runs D6 and D9). 

In this paper we apply the multi mass refinement to 
the cluster 'D' from DMS04. This cluster is well relaxed and 
isolated at z=0 and has an average density profile and inner 
slope close to the mean value. First we mark all particles 
within one percent of the virial radius in the final halo and 
trace them back to the initial conditions. Then we add all 
particles within one comoving Mpc of a marked particle to 
the set of marked particles, and finally we add all particles 
which lie on intersections of any two already marked par- 
ticles on the unperturbed initial grid positions. After these 
two steps there is region with a fairly regular triaxial bound- 
ary which contains only marked particles. The number of 
marked particles grows by almost a factor of 8 during these 
additions, but it is still more than a factor of two smaller 
than the number of particles in the final cluster and a fac- 
tor of ten smaller than the original high resolution volume 
used in DMS04. The computational cost with our code and 
parameters is roughly proportional to the number of high 
resolution particles, therefore we gain about a factor of ten 
with this reduction of the high resolution region. Probably 
one can reduce the high resolution volume further and focus 
even more of the computational effort into the innermost 
region, we plan to explore this possibility with future simu- 
lations. 



2.2 Code and parameters 

The simulations have been performed using PKDGRAV, 
written by Joachim Stadel and Thomas Quinn iStadell200ll) 
using the same cosmological and numerical parameters as 
in DMS04 with a few changes given below and in Ta- 
bic The cosmological parameters are (fl m , SI a, o"8, h) = 



3.0 10 8 






16.2 


0.25 


GJ 


1.0 10 6 


1.8 10* 






13.5 


0.2 


fn 


1.8 10 6 


1.8 10 8 


3-6[p] 


3.8 10 10 


13.5 


0.2 


ill 


1.8 10 6 


1.8 10 8 


38. 6 [p] 


3.8 10 10 


13.5 


0.2 


GJ 


1.8 10 6 


5.2 10 7 






9.0 


0.2 


CD 


6.0 10 6 


5.2 10 7 


15 [p] 


1.4 10 9 


9.0 


0.2 


(i, 


6.0 10 6 


2.2 10 7 






6.8 


0.2 


CI) 


1.4 10 7 


2.4 10 6 


9[c] 


3.0 10 8 


3.3 


0.25 


(2l 


1.3 10 s 


2.4 10 6 


9[P] 


3.0 10 8 


3.3 


0.25 


(I) 


1.3 10 8 


3.0 10 5 


6[c] 


3.75 10 7 


1.7 


0.25 


& 


1.0 10 9 



(0.268,0.732,0.7,0.71). The value of <x 8 = 0.9 given in 
DMS04 is not correct: During the completion of this pa- 
per we found that due to a mistake in the normalization our 
initial conditions have less power than intended. This low- 
ers the typical formation redshifts and halo concentrations 
slightly but does not affect the slopes of the inner density 

profiles. 

We use the GRAFICS2 package jBertschingeJl200 ll) to 
generate the i nitial conditions. The particle time-step crite- 
rion Ati < rjyj ejai , where aj is the acceleration of particle 
"i" , gives almost constant time-steps in the inner regions of 
a halo (see Figure 2 in DMS04), but the dynamical times 
decrease all the way down to the center. Therefore the time- 
step criterion was slightly modified, to make sure enough 
time-steps are taken also near the halo centers: Instead of 

Ati < r)^[7la l (1) 
we now use 

At < mm(r t y / e/a i ,ri/4y /r Gp~i) , (2) 

where pi is the density at the position of particle "i", ob- 
tained by smoothing over 64 nearest neighbors. We used 
r\ = 0.25 for runs DM25 and DM50. Note that in the 
inner region of a CDM halo p(r) ~ 0.6p(< r), i.e. 0.8 
y Gp(ri) ~ y/ Gp{< ri) therefore the condition (J2J with 
rj = 0.25 assures that at least 12 time-steps per local dy- 
namical time 1/ 'yj Gp(< r») are taken. 

The time-steps are obtained by dividing the main time- 
step (to/200) by a factor of two until condition @ is fulfilled. 
In runs DM25 and DM50 the smallest particle time-steps 
are to/51200. According to Figure 2 in DMS04 this time- 
step is sufficient to resolve smaller scales than 0.1 percent 
of the virial radius, i.e. less than the limit set by the mass 
resolution, even in run DM50. 

The smaller time-steps in the inner regions of the clus- 
ter are crucial: In Figure Q we compare two runs which 
only differ in the time-step criterion. DM251t was run with 
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Figure 1. Density profiles in physical (not comoving) coordi- 
nates at redshifts 4.4 and 0.8. The two runs have equal mass 
resolution but different time-steps and softening. The arrow indi- 
cates the resolution limit set by the particle mass. The run with 
the larger time-steps and softening underestimates the dark mat- 
ter density outside of the resolution scale. 

the standard criterion and 77 = 0.2, for run DM25 we 
used the more stringent, computationally more expensive 
criterion (|5J and r\ — 0.25. The difference in CPU time is 
about a factor of two. At z=0.8 the densities in run DM251t 
are clearly lower out to 0.003 virial radii which also af- 
fects part of the region we aim to resolve with this run 
(''resolved = 0.0019r v i r ). Due to the high computational cost 
of these runs we cannot perform a complete series of conver- 
gence test at this high resolution but due to the monotonic 
convergence behav ior of PKDGRAV for shorter time-steps 
iPower et al]|2003T) we are confident that DM25 is a bet- 
ter approximation to the true CDM density profile of this 
cluster. 

Our time-stepping test confirms that the time resolu- 
tion in DMS04 was sufficient to resolve the minimum scale 
of 0.3% virial radii set by their mass resolution. For the pur- 
pose of this work, i.e. to resolve a region even closer to the 
center smaller time-steps are necessary. These two runs il- 
lustrate nicely how a numerical parameter or criterion that 
passes convergence tests performed at low or medium reso- 
lution can introduce substantial errors if employed in high 
resolution runs. 



2.3 Testing the multi mass technique 

Reducing the high resolution region in the way described in 
Section ^. ll produces multi mass virialised systems, i.e. halos 
where particles of different mass are mixed up with each 
other. The inner regions are dominated by light particles 
and the region near the virial radius by heavier particles. 
But one will find particles of both species everywhere in the 
final halo and one has to worry if this mixing introduces 
numerical effects, like energy transfer from the outer part 
to the inner part (from the heavy to the light particles) 



due to two body interactions. This could lead to numerical 
flattening of the de nsity profile and make heavy particles 
sink t o the center ijBinnev fc Knebel 120021 : iDiemand et alJ 
l2004al) . 

To check if the multi mass technique works for cosmo- 
logical simulations we re-ran the simulations D6 and D9 
from DMS04 using a reduced high resolution region. We 
call these multi-mass runs "DM6se" , "DM61e" and "DM9" 
(see Table 0. The next heavier particles in the surround- 
ing region are 216 times more massive in DM6se and DM61e 
and 27 times more massive in DM9. The heavier particles in 
DM61e and DM9 have larger softening to suppress discrete- 
ness effects while DM6se uses the same small softening for 
both species. Figure [3] shows that the density profiles of the 
fully refined run D9 and the partially refined run DM9 are 
identical over the entire resolved range. Figure shows that 
the same is true for run DM61e, the larger mass ratio of 216 
does not introduce any deviation form the density profile of 
the fully refined run. 

A small softening in the heavier species (run DM6sl) 
does introduce errors in the final density profile (Figure 0. 
The total mass profile is shallower near the resolved radius 
and has a high density bump below the resolved scale. The 
light particles are more extended and the bump is caused by 
a cold, dense condensation of six heavy particles within 0.004 
r v i r . These six heavy particles have a 3D velocity dispersion 
of only 273 km/s, while the light particles in the same region 
are much hotter, a^o = 926 km/s. They are hotter than the 
particles in the same region in run D6 and DM61e (both 
have only light particles in this inner part), the dispersion 
are 722 km/s for D6 and 708 km/s for DM61e. 

These tests indicate that the reduced refinement regions 
work well in runs D9M and DM61e and therefore we used the 
same refinement regions to set up the higher resolution run 
DM25. In this run the heavier particles are 125 times more 
massive than the high resolution particles and they have a 
softening of 9 kpc. For run DM50 we refined only the inner 
part of the most massive cluster progenitor at z—4.4 in the 
same way as the final cluster in runs DM61e, DM6se, DM9 
and DM25. In run DM50 the heavier particles are also 125 
times more massive than the high resolution particles. 

Figure [3] shows how the initially separated species of 
light and heavy particles mix up during the the runs DM9, 
DM25 and DM50. The density profiles profiles of DM61e 
and DM9 do not suffer from numerical effects due to the 
multi-mass setup. This indicates that the same is true for 
run DM25 which has the same refinement regions. In run 
DM50 the amount and location of mixing at z=4.4 relative 
to r2oo is very similar to the situation if DM9 at z=0.0, 
therefore we expect DM50 to have the same density profile 
as a fully refined cluster, i.e. as a cluster resolved with a 
billion particles. 



3 THE INNER DENSITY PROFILES 

Here we try to answer the question if the inner density 
profiles of dark matter halos have a constant density or 
a cusp p(r) oc r '. At resolutions of up to 25 million 
particles within the virial radius there is no evident con- 
verge nce toward any constant inner slope iFukushiee et all 
120041 DMS04). 
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Figure 3. Tests of multi-mass refinement and convergence. The upper left panel shows that run D9 which contains only high resolution 
particles within the virial radius has the same density profile as the multi-mass run DM9. The z=0.8 profiles are shifted downward by a 
factor of ten for clarity. The arrows indicate the convergence radius of run D9 estimated in DMS04. The lower left panel shows the high 
and low resolution particles in run DM50 at z=4.4. The panels on the right illustrate the mixing of light and heavy particles in runs 
DM9 and DM25 which have the same refinement regions. 



3.1 Results of run DM25 

Run DM25 has an effective resolution corresponding to 127 
million particles within the virial radius and a force resolu- 
tion of 0.48 x 10 _3 r v i r . At this up to now unmatched reso- 
lution the inner slope is roughly constant from the resolved 
radius (see Figure |1J out to about one percent of the virial 
radius of the final cluster. 

Run D12 resolves the same cluster with 14 million par- 
ticles and shows no convergence to a constant inner slope. 
Note that the "D" cluster is one of six clusters analyzed in 
DMS04 and its inner profile is not special and rather close 
to the sample average. 

Figure |1] indicates that there is a cusp in the centers of 
cold dark matter clusters and it becomes apparent only at 
this very high numerical resolution. The non-constant slopes 
just near the convergence scale are probably due to the first 
signs of numerical flattening that set in at this scale. At 
higher densities below the resolved scales one cannot make 



any robust predictions yet, but if one has to extrapolate into 
this region Figure 2] motivates the choice of a cusp p(r) oc 
r~ 7 with 7 ~ 1.2. 

3.2 Resolving the very inner density profile at 
z=4.4 (run DM50) 

Mass accretion histories show that the inner part of 
CD M halos is a ssembled in an early phase of fast accre- 
tion Jvan den Boschll2002l : IWechsler et~al"ll2002t IZhao et alJ 
2003) and recent high resolution simulations revealed that 
the inner density profile does not evolve at low re dshift 
jFukushiee Kawai fc Makinol Eool iTasitsiomi et alJ 120041: 
iReed et al.ll2005D . Figure 0] confirms that the inner density 
profile of runs D12 and D5 do not change from z=0.8 to 
z=0. 

Therefore in run DM50 we focus our computational ef- 
fort even more on the early evolution of the inner profile. 
We refine the inner region of the most massive progenitor 
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Figure 2. Tests of the multi-mass refinement technique. The 
upper three lines shows the total density profile at z=0 from the 
fully refined run D6 (solid lines) and the multi-mass runs DM6se 
(dashed) and DM61e (dashed dotted). The lower lines (same line 
styles, offset by two magnitudes for clarity) show the density pro- 
files of the two particle species, i.e. of the light ones (lines with- 
out symbols) and of the heavier ones (lines with symbols: filled 
squares for D6se, open circles for D61e). The vertical dashed line 
indicates the innermost resolved scale. In the multi-mass run with 
more softened heavier particles (D61h) the inner profile is dom- 
inated by light particles and identical to the fully refined run 
of the same cluster (D6). When the heavier particles have short 
softenings some of them spiral into the center due to dynamical 
friction and transfer heat to the light particles. This affects the 
total density profile, i.e. it is lower near the resolved scale and 
has a bump due to a condensation of cold, massive particles very 
close to the center. 



identified in run DM25 at z—4.4. Since the refinement region 
needed is much smaller than the one of DM9 or DM25 and 
we only run the simulation to z=4.4 it is feasible to go to a 
much better mass and force resolution. The high resolution 
particles in run DM50 are a billion times lighter than the 
final cluster. 

Figure [£] shows that the density profile of run DM50 
at z=4.4 is cuspy down to the resolved radius (0.1 % of 
the final virial radius). As in run DM25 the slopes begin to 
shallow just at the converged scale due to numerical flat- 
tening. The profile of DM50 at z=4.4 supports the finding 
from run DM25 that the inner profile follows a steep power 
law p tx r~ 1,2 . At the higher resolution of run DM50 we find 
substantially higher physical densities in the cluster center 
at z=4.4 compared to lower resolution runs like DM25. This 
suggests that a run like DM50 evolved to low redshift would 
also yield substantially higher central densities as currently 
resolved in the centers of runs like D12 and DM25. 



3.2.1 Estimating the z=0 profile of a billion particle halo 

Now we go one step further and use the information from all 
the "D" -series runs to try to estimate the density profile one 
would obtain if one simulates this cluster with a billion par- 
ticle all the way to present time, a run which would be pos- 



Figure 4. Logarithmic slope of the density profile of run DM25 
at z=0.8. The slopes of runs D5 and D12 at z=0.8 and z=0 are also 
shown for comparison. The arrows indicates the estimated con- 
vergence radii. Note that although the densities at the converged 
scales are within 10 percent the density gradients can already be 
substantially smaller. 




Figure 5. Logarithmic slope of the density profile of run D5, 
DM25 and DM50 at z=4.4. The arrows indicates the estimated 
convergence radii. A constant inner slope of about -1.2 is evident 
in the highest resolution run DM50. The increase of the slopes 
around the resolved radii is due to the onset of numerical flatten- 
ing. 



sible but extremely expensive with today's computational 
resources. From Figure |S| one finds that the density profile 
of run DM25 near its resolution scale shifts upward by a con- 
stant factor of 1.4 from z=4.4 to z=0.8. The density around 
0.01 r v i r ,z=o is constant form z=0.8 to z=0, see run D5 in 
Figure The inner density profile slopes are constant even 
longer, i.e. from z=4.4 to z=0, see Figures0|and|S| Therefore 
we estimate the z=0 profile of run DM50 by rescaling the 
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Figure 6. Density profiles in physical (not comoving) coordi- 
nates at redshifts 4.4, 0.8 and 0. Arrows mark the resolved scales 
of each run. The densities in the inner part do not evolve between 
z=0.8 and z=0 and the inner slopes remain constant even from 
z=4.4 to z=0. Using these observations we are able to estimate 
the final profile of a billion particle halo (upper solid line). 



z=4.4 profile of DM50 by a factor 1.4 and using the z=0 pro- 
file of run D12 outside of 0.005 r v ir,z=o (see Figure^J. The 
extrapolated z=0 profile of run DM50 should be regarded 
as a best guess for the density profile of an average CDM 
cluster resolved with a billion particles. A (multi-mass) sim- 
ulation with this (effective) resolution evolved to redshift 
zero would be needed to check the accuracy of the estimate 
performed here. Note that our conclusions are based on the 
z = 0.8 results from run DM25 and not on the somewhat 
uncertain z — extrapolation proposed in this section (but 
they are consistent with it). 



3.3 Inner slope estimates based on the enclosed 
mass 

For a mass distribution which follows p(r) oc r -7 all the 
way in to r = the slope 7 can be calculated at any ra- 
dius using the local density and the mean enclosed density 
jNavarro et al.ll2004h : j*(r) = 3(1 -p(r)/p(< r)). For simu- 
lated CDM density profiles where 7 becomes smaller towards 
the center 7*(r) is an upper limit for the asymptotic inner 
slope as long as both p(r) and p(< r) have fully converged at 
radius r. Convergence tests show that the enclosed density 
p(< r) converges slower than the local density and p(< r) is 
generally under estimated near r res olved due to missing mass 
within r rcso i vcd JPower et alJl2003l: DMS04). Figure |7| shows 
7*(r) for the two highest resolution runs available at z=4.4 
and z=0.8. We also plot the fractions of the local densities of 
the two runs and the fractions of enclosed densities to illus- 
trate the different convergence scales of local and cumulative 
quantities. Figure [7| confirms that at the estimated resolved 
scales for D12 and DM25 the local densities are within 10% 




DM25 



0.0005 0.001 0.005 0.01 

Figure 7. 7* (r) for the two highest resolution runs at z=4.4 and 
z=0.8 (solid lines) and fractions of the densities of these two runs 
(dotted lines for p(r) and dashed lines for p(< r)). Due to dif- 
ferent convergence rates in local and cumulative quantities 7* (r) 
values from the lower resolution runs lie below the higher resolu- 
tion results in the inner part of the halo. The arrows at r reso i ve( j 
correct for this effect based on the following observations: The 
ratio p(r reso i ve( j)/p(< ^resolved) i s typically underestimated (0.87 
of the high resolution value) due to a small deficit in local density 
(0.95 of the true value) and a larger one in the enclosed density 
(0.83 of the true value) due to missing mass in the innermost re- 
gions. Underestimating p(r reBolved )/ p(< r rcso i vcd ) by 0.87 leads 
to 7* values which are to small by about 0.3. 

of the higher resolution runs 2 . The typical differences are 
even smaller (about 5%). The enclosed density p(< r) how- 
ever converges slower: At JVesolved we find that the values 
are only about 0.83 of those measured in the higher reso- 
lution runs. This causes the ratio p(r re3 oived)/p(< n-esdved) 
to be underestimated (about 0.87 of the true value). This 
propagates into a larger relative error in 7* (r reso ived) which 
turns out to be too low by about 0.3 for the profiles studied 
here (given the arrows at r re soived in Figure |7J. The different 
convergence rates of local and cumulative quantities tend to 
produce artificially low 7* (r) values and this effect becomes 
especially large near r rcso i vc d. The significance of 7*(r) ap- 
pears to be difficult to interpret, but the convergence tests 
presented here and in DMS04 suggest that 7* (r rcso i vc( j) is 
not a robust upper limit for the asymptotic inner slope. 

3.4 Cored and cuspy fitting functions 

In this section we fit one cuspy and two recently proposed 
cored functions to the density profiles of DM25 at z=0.8 and 
to the tentative z=0 extraploation from run DM50. From the 
last section we expect the cuspy function to work better in 

2 We determine r reS olved by demanding that the local density has 
to be within ten percent of the value from a much higher resolu- 
tion run and in cases where no such run is available the measured 

convergence radii form lower resolution runs are rescalcd using the 

— 1/3 

mean inter-particle separation resolved oc N . (see DMS04). 
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Figure 8. Density profile of run DM25 at z=0.8 and fits with 
three different functions. 
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Figure 9. Density profile of run DM50 extrapolated to z=0.0 
and fits with three different functions. 



the inner part but we try to fit also the cored profiles for 
comparison. 

We use a general Q/37-profile that asymptotes to a cen- 
tral cusp p(r) oc r ': 



Ps 



(r/r s )T(l + (r/r a )-)05-^)/ c 



(3) 



If one takes a, f3 and 7 as free parameter one encounters 
strong degeneracies, i.e. very different combinations of pa- 
rameter values can fit a typical density profile equally well 
jKlvpin et al]|200lD . Therefore we fix the outer slope /3 = 3 
and the turnover parameter a = 1. For comparison the 
NFW profile has (a, j3, 7) = (1,3,1), the M99 profile has 
(a,f3,y) = (1.5,3,1.5). We fit the three parameters 7, r s 
and p„ to the data . 

iNavarro et alJ i2004Tl proposed a different fitting func- 
tion which curves smoothly over to a constant density at 
small radii: 



HMr)/ps) = (-2/a N ) \{r/r s ) c 



1] 



(4) 



qn determines how fas t this profile tu r ns aw ay from a power 
law in the inner part. INavarro et alJ a2004Tl found that qn 
is independent of halo mass and qn = 0.172 ± 0.032 for all 
their simulations, including galaxy and dwarf halos. 

Another profile that also curves away fro m power law 
behav ior in the inner part was proposed by IStoehr et alJ 
(2002): 



PSWTS(0 



Vmax jq-2«SWTS [ l0g ( Tma as ) ] _L X 

4nG r 2 



x 



4ti loi 



\ 1 max / 



where Vmax is the peak value of the circular velocity, r max 
is the radius of the peak and oswts determines how fast 
the profile tu rns away from an power law near the center. 
IStoehrl i2004T) found that cluster profiles are well fitted with 
this formula using oswts values between 0.093 and 0.15. 

These three functions were fitted to the data from z=0.8 
by minimizing the relative density differences in each of 



about 20 logarithmically spaced bins in the range resolved 
by DM25 (i.e. form 0.0019r v i rjZ =o = 3.3 kpc to r v i r ,z=o = 
1750 kpc). At z=0 we use the resolved range of D12 for 
the fits (i.e. form 0.0039r v j r ,z=o = 6.8 kpc to r v i r ,z=o). The 
resulting best fit values and the root mean squares of the 
relative density differences are given in Table |21 

At z=0.8 the average residuals of the three fits are sim- 
ilar, but they are dominated by the contribution from the 
outer parts of the cluster (see Figure 6 in DMS04). Figures 
|H] and |5] show that in the inner part the cuspy profile de- 
scribes the data better. Both cored profiles underestimate 
the measured density at the resolution limit both at z=0.8 
and in the estimated z=0 profile. These profiles lie below the 
measured density profiles even inside of r rea oivcd where one 
has to expect that the next generation of simulations will be 
able to resolve even higher densities. 

Fi eures [TU1 and [Til show the slopes of the simulated pro- 
file in comparison with the slopes of the best fits. Again it 
is evident that in the inner part the cuspy profile describes 
the real density run better. 



4 CONCLUSIONS 

The main conclusions of this work are the following: 

• It is possible to use different mass particles to resolve 
one halo in cosmological CDM simulations without affecting 
the resulting density profiles. 

• This "multi-mass" technique allows a reduction of the 
necessary number of particles and the computational cost by 
at least one order of magnitude without loss of resolution in 
the central region of the halo. 

• We confirm that the inner profile of a typical CDM 
cluster does not evolve since about redshift one. 

• The logarithmic slope of the dark matter density pro- 
file converges to a roughly constant value in the inner part 
of cluster halos. This probably holds also for smaller sys- 
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Table 2. Density profile parameters of run DM25 at z=0.8 and of DM50 extrapolated to z=0. A is the root 
mean square of (p — pat)/p for the three fitting functions used. 



redshift 


7G 


r sG [kpc] 


A G 


a N 


r s N [kpc] 


A N 


a SWTS 


fmai SWTS [ k P c ] 


Aswts 


0.8 


1.20 


260 


0.075 


0.157 


233 


0.076 


0.130 


565 


0.087 


0.0 


1.20 


283 


0.059 


0.162 


236 


0.133 


0.140 


518 


0.179 



• The cluster studied here has a central cusp p cx r~~' with 
a slope of about 7 = 1.2. From earlier studies (DMS04) we 
expect this inner profile to be close to the average and the 
scatter is about 0.15. 

• P rofiles with a core iStoehr et alJ l20ol3lNavarro et alJ 
120041) underestimate the measured dark matter density at 
(and even inside of) the current resolution limit. 



Data 

Navarro et al 03 

- a|iy-profile with 

best fit inner slope y = 1 .2 
(a=1 , |i=3) 

SWTS 
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Figure 10. Logarithmic slopes of the measured and fitted den- 
sity profiles from Figure 151 
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Figure 11. Logarithmic slope of the extrapolated z=0 DM50 
density profile and of the fitted density profiles from Figure l9l 

tems (like galaxy and dwarf halos) but there it is even more 
difficult to numerically resolve the cusps. 

• At resolutions around 10 million particles per halo the 
inner slope appears to approach zero continuously but this 
impression is caused by numerical flattening of the profiles 
due to insufficient mass resolution. 
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